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Abstract. This document is a companion for the Maple program 
Summing a polynomial function over integral points of a 
polygon. It contains two parts. First, we see what this programs 
does. In the second part, we briefly recall the mathematical back- 
ground. 



1. Introduction 

The present article is a user's guide for the Maple program Sum- 
ming a polynomial function over integral points of a polygon, 
available at http : //www. math. polytechnique . f r/~berline/maple .html. 

The Maple program contains two types of computation. The first com- 
putation does just what the title says. The input consists of a finite 
set of rational points in Q^, whose convex hull is a polygon p, and a 
polynomial h{x, y) with rational coefficients. The output is the sum 

hix,y). 

The second computation returns the function of t G N which arises 
when the polytope p is dilated by t. 

E{t):= Yl ^(^'?/)- 

{x,y)&pnZ^ 

This function is a quasi-polynomial, meaning that is has the form 

dcg h+2 

E{t)= J2 

where the coefficients depend only on t mod q, where q is the smallest 
integer such that gp has integral vertices. The function E{t) is called 
the weighted Ehrhart quasi-polynomial of p with respect to the weight 
h{x,y). 
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We apply two methods, the first one for a fixed polygon, the second 
one for the computation of the weighted Ehrhart quasi-polynomial. 
The first method is based directly on Brion's formula ([2]), [3], while the 
second method is based on the local Euler-Maclaurin formula of |2j. 
Both methods use Barvinok's decomposition into unimodular cones 
[1]. Although they are very similar, the first method is faster when we 
deal with a fixed polygon, while the second is faster when we want the 
Ehrhart quasi-polynomial. 

The software libraries LattE [4j (improved version in [5j) and Barvi- 
nok include the computation of the number of points of a ratio- 
nal polytope in any dimension, together with many other applications. 
Moreover, the weighted Ehrhart polynomials in any dimension are com- 
puted in Barvinok. The present program, in dimension two, is based 
on the same principles: Brion's formula and Barvinok's decomposition 
of cones. We use however some new ideas on "renormalisation" of Lau- 
rent series from [2] to speed up the computation. In the future, we will 
generalize it to higher dimensions. 



2.1. Summing a polynomial function over the set of integral 
points of a polygon. Let P C be a finite set of points. Let p C 
be the polygon obtained as the convex hull of the set P. The program 
computes the sum 



of the values of a polynomial h{x,y) over the set of integral points 
contained in p. In particular, when /i = 1, it computes the number of 
integral points in p. 

For a single monomial h{x, y) = x"^^y"^'^, the command is 

>sum_mononiial_polygon(P,ni) ; 

Here P is a set of pairs of rational numbers, and m = [mi, 7712] is a pair 
of non negative integers, the multidegree of the monomial 

If we want just the number of integral of integral points, we can use 
the command 

>number_points_polygon(P) ; 

This number can be also obtained by the command 

>suin_mononiial_polygon (polygon, [0,0] ) ; 

We compute the sum of a polynomial h{x, y) by the command 

>sum_polynomial_polygon(P,h) ; 



2. Main commands 





(x,y)Gpnz2 
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Here P is a set of pairs of rational numbers, and h = hmx"^^y"^^ is 
entered as an expression in x,y. 

Example 1. P is the square {[0, 0], [1, 0], [1, 1], [0, 1]}. 
>square:= {[0,0] , [1,0] , [1,1] , [0,1]}; 

>number_points_polygon(square) ; 

4 

The sum of values x^y^ over the 4 integral points in the square is 

>suin_monomial_polygon (square, [5,5 ] ) ; 

1 

Example 2. Here P is a randomly chosen set of 15 points. 

>P := {[77/8,97/59] , [93/44,70/29], [0,25/12], [25/32,29/48], 
[92/41,57/91], [9/4,1/7], [64/43,31/75], [91/17,33/86], [12/37,77/8], 
[8/5,41/27], [80/67,11/9], [16/73,11/89], [41/20,43/88], 

[32/49 , 59/23] , [77/94 , 65/46] } 

The number of integral points in the convex hull is 45. 
>number_points_polygon(P) ; 

45 

The vertices of the convex hull p of P (listed in counter-clockwise order) 
are obtained with the command: 

>vertices_in_counter_clock_order : =proc (polygon) 

>vertices_in_counter_clock_order (P) ; 
[[0,25/12], [16/73,11/89], [9/4,1/7], [91/17,33/86], [77/8,97/59], 

[12/37,77/8]] 

We compute the sum of x^^y^^ over the set of integral points (x, y) of 
the convex hull p of P. 

>suin_monomial_polygon(P, [32,32] ) ; 

987532646688766560932727042325214847653263886 

We compute the sum of x^'^y^'^ + 7 over all the integral points (x, y) of 
the polygon p. (the preceding number +7 times 45) 

>h:= x~{32}y~{32>+7; 
>suin_polynoniial_polygon(P,h) ; 

987532646688766560932727042325214847653264201 

2.2. Weighted Ehrhart polynomial of a polygon. Our program 
computes also the weighted Ehrhart quasi-polynomial of a polygon. 
For brevity, we treat only the case where the weight is a monomial 
h{x,y) — x^^y'^'^. When the polygon is dilated by a non negative 
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integer t, and if g is a positive integer such that qp has integral vertices, 
the function of t given by 

t ^ J2 ^"''y"'' 

is a quasi-polynomial S{t) = Yl^o'"^^^'^ Ei{t)f of degree mi + m2 + 2. 
The coefficients Ei{t) are functions of t modulo q. This program 
computes these coefficients Ei{t) in terms of the symbolic function 
fmod{p * t, q) which stands for (t i— pt mod q). We can either obtain 
each individual coefficient Ei(t) or the full weighted Ehrhart polyno- 
mial S{t). 

Here are the commands: 

> coeff_t_Ehrhart_polygon(i,t,P,in) ; 

The input consists of i an integer, t a letter , P a set of points and 
m = [mi, 7712] a pair of integers which represents the weight; the output 
is the coefficient Ei{t). 

>Ehrhart_polynoinial_polygon (t , P , m) ; 

Input is as in the previous command, except i is not needed. The 
output is the full Ehrhart polynomial S (t) . 

Examples 

>transsquare : ={ [-1/2 , -1/2] { [1/2 , -1/2] , [1/2 , 1/2] , [-1/2 , 1/2] } ; 

> coeff_t_Ehrhart .polygon (0,t, square, [0,0]); 

1 

> coeff _t_Ehrhart_polygon(0,t,transsquare, [0,0]); 
-2*fmod(t, 2)+3/2+l/2*fmod(t, 2)~2 

> Ehrhart _polynomial_polygon(t , square , [0,0]); 

1 + 2 t + t^2 

> Ehrhart _polynomial_polygon(t,transsquare, [0,0]) ; 

(fmod(t, 2)-l)"2+(-2*fmod(t, 2)+2)*t+t"2 

2.3. Experiments. The following experiments were done with a lap- 
top, processor 1,86 GHz, RAM 782 MHz,0,99 Go. 

>A:={ [ (-567337) /102495, -1414975/95662] , [1/3,1/5] , [-88141/20499,12732/47831]}; 

>largeA: ={[1000* (-567337) /102495, 1000* (-1414975/95662)] , 
[1000*1/3,1000*1/5] , [-1000*88141/20499,1000*12732/47831]}; 

> number_points_polygon(A) ; 

36 
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> niunber_points_polygon(largeA) ; 

34922612 

In the next experiments, we indicate the time of computation T in 
seconds. The number of integral points in the rational triangle with 
vertices A is 36. If we dilate A by the factor 1000, we obtain the triangle 
largeA where the number of points (34922612) is approximatively 10^ 
times larger. Observe that we compute in 14 seconds the sum of the 
large degree monomial y) = x^^y^"^ over the set of integral points of 
A, and that we compute in 16 seconds the sum of the same monomial 
over the integral points of largeA. The computation time is almost the 
same, although any computation by enumeration would be 10^ times 
longer. 

> T:=time(): suin_monomial_polygon(A, [32,32] ) ;Time:=time()-T; 
11156693714080121436809683716369682546812787494001398139657 

Time := 1.766 

> T : =t ime ( ) : suin_monoinial_polygon (A , [64 , 64] ) ; Time : =t ime ( ) -T ; 
10691662746975383171690687952963005219723639375189814 
217756070191566530558879\ 

3836513555847334896253718879462978590217 

Time := 13.640 

> T:=time(): sum_monomial_polygon (largeA, [64,64] ) : Time :=t ime ()-T; 
17831035913722066043589677840496661987989193563450057671832979767102708226068\ 

905195223428957659882216123374803724362290728944933635792703052976782671238\ 
401601191375977184037799597789861617132380131198911864015293592136365221852\ 
952449214916133197928922419462989960495593699297670700652853834584439172901\ 
857916119694620105996329573478014513449383738873972550889051937620201341771\ 
829110756841837358870588454172079624770000592845281113102517836579429050870\ 
20099703621578931359063825440122383120351301766010118556183 

Time:= 15.516 



Finally, we computed the weighted Ehrhart polynomial with weight 
^,32^32 Q^gj. ^jjg triangle with vertices[[-567337/102495, -1414975/95662], 
[88141, 292844676/6833], [-88141/20499, 12732/47831]]. We compute 

the coefficient of for example. The time of computation is 268 sec- 
onds. The result is too big to be printed here, as it involves many 
functions fmod{c*t, D) where D runs though the denominators of the 
coordinates of the vertices of A (large numbers) . 

> T:=time() :coeff_t_Ehrhart_polygon(2,t, [[-567337/102495, 
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-1414975/95662], [88141,292844676/6833], [-88141/20499, 
12732/47831]] , [32,32]) : Time : =time () -T; 

Time := 268.266 

3. Mathematical background 

The first method is for a fixed polygon, the second one for the com- 
putation of the weighted Ehrhart quasi-polynomial. 

3.1. First method: Brion's formula, Barvinok's decomposition 
into unimodular cones and iterated Laurent series. Let p be a 

convex polygon in with rational vertices Si,l < i < n + 1. We want 
to compute the sum 



(1) Yl ^'"'^ 



We start by observing that ([T]) is equal to the coefficient of ^ , , in 



Our method is based on Brion's formula ([2]). Brion's formula is the 
generalization of the following formula for the sum of geometric pro- 
gressions over the interval [A, 5] (with A< B integers): 

^ pBi 

A 

For any rational polygon q C define 

This a meromorphic function near ^ = 0. Moreover the map q ^ 
S'(q)(^) is a valuation on the set of rational polyhedra, and 5'(q) = if 
q contains a line. Brion's formula is the following. Let Cj be the cone 
at vertex Si of the polygon. 

71+1 

(2) S{^) = Y,S{U). 

i=l 

Each term S{ti){^) = J^xeciCiZ^ ^^^'^^ in ([2]) is a meromorphic function 
near ^ = 0. The poles cancel and the sum is a holomorphic function 

of ^. Thus we compute ([1]) as the coefficient of ^^|^_^, in the right- 
hand-side of ([2]). We actually compute the individual contribution of 
each cone Cj (associated to the vertex Si) to the sum. The coefficient of 
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in the meromorphic function S'(Cj) of two variables .^i, ^2 has no 
intrinsic meaning. Our method consists in applying iterated Laurent 
series expansions to S'(Cj)(^) with respect to the variables ^1 then ^2- 
We obtain a Laurent series L{ci) in the ring Q[^i, ^2, ^2""*^] we 
compute the coefficient ^ , . in L(cj. 

Thus, in order to compute the contribution of a vertex s to the sum 
(I2I), we need to compute S'(c)(^) for the supporting cone c. The crucial 
tool here is Barvinok's decomposition into unimodular cones. Actually, 
we use the following variant of Barvinok's decomposition, (procedure 
signed_decomp). 

Let c be a simplicial cone in M'^. Let V^,, for i = 1,. . . ,d, be the 
generators of c. Let K be a vector in M.'^. We write V = Ylii'^i^i- We 
split [Vi, . . . , Vrf] into three parts, as follows. 

L+ := [Xi, . . . , Xfc] 
formed by the Vi such that > 0, 

L„ := [Yi, . . . , 1^] 

formed by the Vi such that < 0, 

'■= {Zi, . . . , Zb} 

formed by the Vi such that Ui = 0. 

Then we have the equality of characteristic functions modulo char- 
acteristic functions of cones containing lines. 

k 

(-l)(^+i)[c] = ^(-l)^+i[c(Xi, . . . , X,_i, -X,+i, . . . , -X,, V, L_, Lo)] + 

i=l 

m 

5^(-iy+'=[c(L+, -V, -Fi, . . . , . . . , r^, Lo)]. 

Remark. This decomposition is not the stellar decomposition. It 
involves only cones of maximal dimension d. It avoids the dualizing 
trick of Br ion. 

Example, c = M+ei ©M+e2, V = ei + 62, so that L_ and Lq are empty 
and k = 2. Then 

-[c] = c{V, -62) - c(ei, V) - c[e2, -62, ei]. 

Indeed [c(e2, —62, ei)] — [c] is equal to the characteristic function of the 
quadrant (ci, —62) minus that of the half- line M+Ci. This is also the 
case for [c(V, -62)] - [c(ei, V)]. 
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If we use a lattice vector V with sufficiently small coordinates in 
the basis (T^), the cones appearing in this decomposition have indices 
smaller than c. One obtains such a short vector V by the Lenstra- 
Lenstra-Lovasz algorithm. By a repeated application of this decompo- 
sition, one obtains a decomposition of c in a signed sum of unimodular 
cones Cz (modulo cones containing hues). As S{a) = for a cone a 
which contains a line, we can use this decomposition to compute S{c). 

For a unimodular cone c, the sum S{c) has a simple closed expression. 
Let (Vi, V2) be primitive generators of the edges of c and let s be its 
vertex. Let s be the unique integral point contained in the semi-closed 
box 

{s + tiVi+t2V2,0 <ti < 1} 

If s = si^i + S2V2, then s — siVi + S2V2 with Sj = ceil{si). Then 

(3) '^('^)(0= (i_g(e,yo)(i_e(«,V2>)- 

In order to simplify the computation of iterated Laurent series, we 
introduce the analytic function 

n=0 ^ ' 

where b{n, u) are the Bernoulli polynomials. Writing 

uX Y 

(4) —-^^B{X,u)- 



we obtain 
where 



%)(0=^ + G + it:, 



B{{i,v,)rsi)B{{i,v2)rs2) 

is an analytic function of ^, 

1 

R :-- 



(^,^l)(e,^2>' 

We replace j^^^ and by their iterated Laurent series expansion 
in the ring R[^i, ^, ^2, ^2^^]- For example, if Vi — [2, 1], we write 



1 



26 + 6 6(1 + 26/6) 6,^0 



-^(-i)'=2'=(6/6) 
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We then replace S'(c) by the corresponding element in Q[[^i, ^2, ^ ^]] 
and we take the coefficient of ^"[^^CT^. 

Remark The weighted Ehrhart polynomial can also be computed by 
this method. We did not write the corresponding algorithm in the 
Maple file, because we observed that a faster algorithm is given by the 
second method which we describe in the next section. However, let us 
explain what one should do. When the polytope p is dilated in tp, its 
vertices are dilated by t, while the edges of the cones at vertices do not 
change. Thus we have to compute 

= (i_e(«,Vi))(i_e(?,y2>) 
where now St is the unique point with integral coordinates in the box 

{ts + UiVi + U2V2, 0<Ui< 1}. 
If s = siVi + S2V2, with Si = Pi/qi, we see that 

St = [tsi + mod{-tpi,qi)/qi,ts2 - mod{-tp2, q2) / q2]- 

The iterated Laurent series in ^1,^2 has coefficients which are poly- 
nomials in t and the periodic functions mod{tpi, qi) . We extract the 
coefficient of t^C'C- 

3.2. Second method. Weighted Ehrhart quasi-polynomial us- 
ing local Euler-Maclaurin formula. We now recall the results of 
[2] and explain how they can be applied to the computation of the 
weighted Ehrhart quasi-polynomials. Let p be a convex polytope in 
M"^, with rational vertices. Let h{x) be a polynomial function of degree 
r on M"^. We want to compute the sum ^a-gpn^ ^(•^) values h{x) over 
the set of integral points of the polytope p. 

The local Euler-Maclaurin formula has the following form. 

(6) M^)= E [D{P,f)-h 

where ^(p) is the set of all faces of p. For each face f, D(p, f) is a differ- 
ential operator of infinite degree with constant coefficients associated 
to f. The operator -D(p, f) is local, in the sense that it depends only on 
the intersection of p with a neighborhood of any generic point of f . The 
integral on the face f is taken with respect to the Lebesgue measure on 
< f > defined by the lattice Z'^ fl lin(f). Here < f > is the affine span 
of the face f and lin(f) is the linear subspace parallel to < f >. 

Let us recall the construction of the operators D{p,f). We denote 
by t(p, f) the transverse cone to p along f. Using the standard scalar 
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product, t(p,f) is described as the following affine cone in M*^. Let 
lin(f)-'- be the vector subspace orthogonal to lin(f). Then t(p,f) is the 
orthogonal projection on lin(f)-'- of the supporting cone of p along f. 
The operator -D(p, f) is defined in terms of the transverse cone t(p, f), 
as follows. 

For every rational affine cone a C V^, we construct in [2J an analytic 
function ^ /i(a)(^) on R'^. This construction depends on the choice 
of a scalar product. Here we use the standard scalar product. These 
functions /i(a) have nice properties which play a crucial role in our 
method. First, the assigment a ^— /i(a) is a valuation on the set of 
affine cones with a given vertex. Second, it is invariant under lattice 
translations. Furthermore, /i(a) = if a contains a line. 

We define 

D(p,f) = D(Mt(p,f))) 
as the differential operator of infinite degree with constant coefficients, 
with symbol /u(t(p, f))(0- other words, if ^ = (^i, . . . , ^rf), we obtain 
D{p,f) by replacing by ^ in the Taylor series of /i(t(p, f))(0- 

For any positive integer t, we consider the dilated polytope tp and 
the corresponding sum 

S{tp,h)= J2 h{x). 

xetpnA 

From ([6]), it follows easily that the function t \—>- S{tp,h) is given by 
a quasi-polynomial: there exist periodic functions t ^— -^^(p, h,t) on N 
such that 

d+r 

(7) S{tp,h) = Y,Ei{P,h,t)t' 

i=0 

whenever t is a positive integer. Moreover the coefficients Ei{p, h, t) are 
computed using the functions /i(t(tp, tf)). Indeed, let s be the vertex 
of t(p, f) so that t(p,f) = s + to. Then the dilated transverse cone is 
t(tp, tf) = ts + {q. As a ^— > /i(a) is invariant under lattice translations, 
we have 

Mt((t + g)p,tf)) =/i(t(tp,tf)), 

if q is an integer such that qs is a lattice point for the projected lattice, 
or equivalently, such that g < f > contains a lattice point. Thus, the 
coefficients -E'i(p, t) depend only on t mod g, where q is the smallest 
integer such that qp has integral vertices. 

When a is a unimodular affine cone of dimension 1 or 2, the functions 
/i(a) have an explicit form, in terms of the functions -B(X, u) introduced 
in dlD. 
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Let be a one dimensional affine cone of the form (s + M+)y where 
\^ is a primitive vector and s G Q. We have 

(8) f^im) = B{{^,V),ceUis)-s). 

Let o be a two dimensional unimodular affine cone. Let Vi,V2 be 
primitive generators of its edges, such that det(Vi,V2) = 1. For ^ = 
(^1,^2) £ IR^, let Hi = {^,Vi) , for i = 1,2, be the coordinates of 
^ relative to the dual basis (V^i*,^*). We write the vertex of a as 
SiVi + S2V2 with Si E Q. Let ej = ceil{si) — Si, and let Ci = for 
i = 1,2. With these notations, we have 



(9) /^(a)(0 



gei?/l+e2?/2 



[1 - ey^){i - ey^] 



111 

—B{y2 - Ciyi, £2) H B{yi - (72^/2, ei^ 



yi y2 yiy2 



The function /i(a)(^) is actually analytic, although this is not obvious 
on ([H]). In order to compute the contribution of a vertex s of p to the 
sum ([2]), we need to compute /i(c)(^) when c is the two-dimensional 
supporting cone at s. The crucial tool here is Barvinok's decomposition 
into unimodular cones. The valuation property of a t— > fi{a) makes it 
possible to reduce the computation to the unimodular case, and use ([9]). 
Notice that Q returns a function of the relative coordinates (1/1,1/2), 
which we must convert back to a function of the standard coordinates 
(^1,^2); in order to add the contributions of the various unimodular 
cones in Barvinok's decomposition. Actually, since /i(a) = if the 
cone a contains a line, we use the variant of Barvinok's decomposition 
described in the first method. 
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